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Abstract 

In this paper, we present results obtained from our recent studies on the 
location of the innermost stable circular orbit (ISCO) for binary neutron stars 
(BNSs) in several levels of post Newtonian (PN) approximations. We reach the 
following conclusion at present: (1) even in the Newtonian case, there exists the 
ISCO for binary of sufficiently stiff equation of state (EOS). If the mass and the 
radius of each star are fixed, the angular velocity at the ISCO i?isco is larger for 
softer EOS: (2) when we include the first PN correction, there appear roughly 
two kinds of effects. One is the effect to the self-gravity of each star of binary 
and the other is to the gravity acting between two stars. Due to the former one, 
each star of binary becomes compact and the tidal effect is less effective. As 
a result, Q\sco tends to be increased. On the other hand, the latter one has 
the property to destabilize the binary orbit, and ^?isco tends to be decreased. 
If we take into account both effects, however, the former effect is stronger than 
the latter one, and J?isco becomes large with increase of the 1PN correction: 
(3) the feature mentioned above is more remarkable for softer EOS if the mass 
and radius are fixed. This is because for softer EOS, each star has the larger 
central density and is susceptible to the GR correction: (4) there has been no 
self consistent calculation including all the 2PN effects and only exist studies in 
which one merely includes the effect of the 2PN gravity acting between two stars. 
In this case, the effect has the property to destabilize the binary orbit, so that 
■^isco is always smaller than that for the Newtonian case. If we include the PN 
effect of the self-gravity to each star, fi\sco will increase. 



§1. Introduction 

BNSs like PSR1913+16 are believed to be formed after the second supernova 
explosion in the evolution of massive binary stars. At the formation, orbits of BNSs 
will be highly eccentric because a large amount of the mass of the system is ejected 
in the supernova explosion. If the semi-minor axis of BNSs is less than about a few 
times of solar radius, they will coalesce within the age of the universe ~ 10 10 yr due 
to the emission of gravitational waves (GWs). In the evolution of the binary the 
orbital radius as well as the orbital eccentricity decreases so that the orbit becomes 
almost circular when the orbital radius is about 10 2 times of the neutron star (NS) 
radius where the frequency of GWs is ~ 10Hz, i.e., in the frequency band sensitive 
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to the kilometer size laser interferometric GW detector such as LIGO, EP VIRGO, B> 
GEO600,!) and TAMA.ffi After that, evolution of BNSs is divided into two phases: 
One is the so called inspiraling phase where the orbital radius of BNSs is sufficiently 
larger than the radius of each NS. In this phase, the hydrodynamic effect of NS is 
not so important for the orbital evolution that each star of binary is approximately 
regarded as a point mass. The emission time scale of gravitational waves (GWs), 
iGW) is still much longer than the orbital period and as a result, the circular orbit 
of BNSs is stable, so that BNSs evolve adiabatically in the inspiraling phase. The 
second phase is the post inspiraling phase in which the orbital radius becomes a few 
times of the NS radius and the hydrodynamic effect as well as general relativistic 
(GR) effects between two stars become very important. Just after BNSs enter this 
phase, i.e., the orbital radius of BNSs becomes smaller than the innermost stable 
circular orbit (ISCO), the circular orbit becomes unstable due to the hydrodynamic 
and GR effects and BNSs come into a final plunging and merging phase. 

GWs in the inspiraling phase have a characteristic feature; the Fourier spectrum 
of the signal behaves as oc /~ 7 / 6 .Ej However, when the BNSs enter the post inspiral- 
ing phase, their wave form will change drastically because the motion of the BNSs 
changes from circular orbit to plunging one. In particular, the ISCO is the charac- 
teristic orbit of such transition n The transition point will be sensitively determined 
by mass- radius relation of NSB>0> and the EOS of NSs may be determined from 
the mass-radius relation. & Hence, the signal around the ISCO will have important 
information on the equation of state (EOS), or in other words, the internal structure 
of NS. This means that if we detect a signal from BNSs at the ISCO, we will be able 
to constrain the EOS by comparing the signal with theoretical templates. For this 
reason, theoretical studies for the ISCO are urgent. 

In Fig.l, we schematically show the location of the ISCO in each level of ap- 
proximation which has been determined from recent studies: In the Newtonian case, 
there does not exist the ISCO if we regard each star of binary as point masses, but it 
appears for BNSS|_Qf sufficiently stiff EOS if we take into account the hydrodynamic 
effect on BNSs. H> M> This is because each star is tidally deformed by the gravity from 
the companion star and the tidal deformation induces the attractive force which over- 
comes the centrifugal force at small radius. For soft EOS, the degree of the central 
condensation is high, so that the effect of tidal force is not so important compared 
with the stiff EOS case. 

When we include the first PN (1PN) correction of general relativity, not only 
the gravity between two stars, but also the self-gravity in each star become stronger 
than those in the Newtonian case. If the former effect is stronger than the latter one, 
the orbital radius of the ISCO (-Risco) becomes larger than that in the Newtonian 
case, because the former PN effect has the property to destabilize circular orbits as 
we will show later (§3 and §4). On the other hand, if the opposite relation holds, 
Risco becomes smaller than that in the Newtonian case because each star of binary 
is forced to be compact due to the PN correction of the self-gravity and the tidal 
effect is less important. Recently, investigating equilibrium sequences of corotating 
binary iShibata and Taniguchi found that the latter effect is stronger than the former 
oneB©©, and -Risco becomes small due to the 1PN correction. As a result, the 
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angular frequency at the ISCO increases. 

Another approach has been jjerformed by Taniguchi and Nakamura (TN)E-f as 
well as Lai and Wiseman (LW)ll^P and Ogawaguchi and Kojima.H^ Note that the 
methods adopted by three groups are essentially similar. In their approaches, they 
include the second PN (2PN) corrections and some parts of the higher relativistic 
corrections to the gravity between two stars phenomenologically, although they do 
not take into account the PN terms to the self-gravity of each star using the New- 
tonian equilibrium. In this case, the PN attractive force between two stars tends to 
destabilize circular orbits and as a result, -Risco becomes larger while the angular 
velocity becomes smaller compared with the point particle cases. 

Considering these recent developments for understanding the mechanism which 
determines the ISCO, we organize this paper as follows. We analyze equilibrium 
sequences of binary stars assuming that tew is sufficiently larger than the orbital 
period, and argue the stability of binary along the sequence because the ISCO is just 
the critical point of stability. In §2, after we briefly mention why the ISCO exists 
even in the Newtonian case, we describe how the location of the ISCO for fluid stars is 
determined paying particular attention to the incompressible star. In the latter part 
of §2, we also argue the effect of the EOS analyzing equilibrium states of corotating 
binary of compressible stars. In §3 and §4, analyzing equilibrium states of corotating 
binaries, we argue how the location of the ISCO changes due to the 1PN correction 
for incompressible and compressible fluid cases, respectively. In the incompressible 
case, it is possible to obtain a fairly accurate analytic solution of corotating binaries 
even in the 1PN case, so that in §3, we present the method to derive it in detail. On 
the other hand, in the compressible case, we cannot obtain the analytic solution, so 
that we perform numerical computations and show the results in §4. In §5, we review 
the method and results by TN in which the 2PN and some parts of the higher PN 
corrections are taken into account only for the gravity acting between two stars. As 
mentioned above, the results they (and LW) obtained seem to apparently disagree 
with those by the full 1PN study. One of the reasons for this difference is the starting 
point of the problem. In the 1PN case, the zeroth order solution is the Newtonian 
equilibrium and then 1PN effect is taken into account as the correction. While in 
treatment of TN, a semi-relativistic equation of motion (EOM) for point particle 
with ISCO is the zeroth order equation and then the finite size effect of the star is 
taken into account. Another reason is the consistency of the treatment. In the 1PN 
case, 1PN effects of the gravity are consistently taken into account to both the self- 
gravity of each star and the gravity between two stars. While in TN, only 2PN and 
the some part of the higher order effects of gravity between two stars are considered 
phenomenologically. To understand the effect of the self-gravity more clearly, in the 
last subsection of §5, we compute 1PN ISCO removing the 1PN self-gravity force in 
§3. From the results, we suggest that if we include the 1PN effect of the self-gravity 
in TN, i?isco will decrease. We briefly review recent semi-GR calculations in §6. 
Section 7 is devoted to summary. 

Throughout this paper, we use G and c as the gravitational constant and the 
speed of light. We define the Newtonian mass M and the multipole moment Iijk— 
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as 

M = J pd 3 x and hjk— = J pxiXjXk • • • d 3 x, (I'l) 

where p denotes the baryon density of star and the integral is performed inside each 
star of binary. We also define the tracefree part of the quadrupole moment as 



j-»7 = ijj — -Sjj ifcfci (i-2) 
6 k=i 

where 5ij denotes the Kroneker's 5. 

In the following sections, we take the axis for the orbital motion of binary as X3, 
and two stars have equal Newtonian mass M and the same structure for simplicity. 
In §2.1, §3 and §5, we put the center of the star ltd at the origin and center of the 
star 2 at (— R, 0, 0) (i.e., the orbital separation is R and x\ axis is the semimajor 
axis; see Fig. 2), and consider the hydrostatic equilibrium for star 1. Orbital rotation 
is assumed to be counterclockwise in all sections. 

§2. ISCO for binary neutron stars in the Newtonian theory 

For binary of two point particles in the Newtonian theory, it is well known that 
there is no ISCO: In this case, the energy conservation equation can be written as0 

m ( dr\ _ . . . . ml 2 GM t m 

= E ~ V ^ y<r) = -^ — — P-i) 

where m, Mt, £, E, r and t denote the reduced mass, total mass of the binary, the 
specific angular momentum, the energy, the orbital separation, and the coordinate 
time, respectively. The circular orbit of its radius ro = I 2 /GMt is always stable up 
to the contact of binary because d 2 V/dr 2 is always positive at r = tq. 

We can also explain the stability of the circular orbit in the following way: The 
energy of the circular orbit is 

2r 1 ; 

This expression means that if E < 0, the circular orbit always exists up to ro = 0. 
This also shows the existence of the stable circular orbit for < ro < 00. Note that 
this conclusion follows from the fact that the centrifugal force is always stronger than 
the gravitational force at small radii (r < ro). 

This is not necessarily the case for NS case if we remind the finite extension of 
the star and take into account the hydrodynamic effect on each NS. Due to the tidal 
field from the companion star, each NS of BNSs will be tidally deformed. In such a 
case, V(r) is approximately written as 

m£ 2 GM t m 3GM t 

v{r) 27^ ll(t) ' (2 ' 3) 



Here, "the center of mass" indicates that defined in the Newtonian order. We should be 
careful to define it in the PN order (see §3). 
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where we assume that the center of each star is located on the x\ axis and is 11 
component of the tracefree part of the quadrupole moment of each star. Here, we only 
take into account the lowest order multipole among all the multipole moments which 
are generated by the tidal deformation. We note that changes approximately as 
0(r~ 3 ) if the star does not have an intrinsic spin. So that the third term in Eq. (|2-3| ) 
behaves as 0(r~ 6 ) in reality. Thus, the effect of the tidal deformation is very small 
for large r, but it cannot be negligible in the case of small r. If the tidal deformation 
is sufficiently large, it may overcome the centrifugal potential for a sufficiently small 
r, and in such a case, the orbit becomes unstable before each star of binary comes 
into contact; i.e., for E < E CT [ t where .Ed-it is a critical value of the energy, no circular 
orbit exists. This is essentially what Lai, Rasio, and Shapiro (LRS) pointed outtf II2P 
already. 

2.1. Incompressible case 

We explain the qualitative fact mentioned above by showing sequences of equi- 
librium states of incompressible binary stars in more details. In the Newtonian case, 
the hydrostatic equation in the rotating frame of the angular velocity f2 isEJ) 
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dui dP d 
^ ~{ 5 dxj dxi ^ dx - 

where Ui, P, and denote the velocity in the rotating frame, the pressure, and the 
completely antisymmetric tensor, respectively. U is the Newtonian potential which 
satisfies 

AU = -4ttG P , (2-5) 

and it can be split into two parts; one is the self-gravity part, and the other is the 
contribution from the companion star. We separate U as 

U = U 1 ^ 1 + U 2 ^ 1 , (2-6) 

where U 1 ^ 1 denotes the self-gravity part and U 2 ^ 1 denotes the contribution from 
the companion. These are written, respectively, as 

[/i-i= ttG^-EW)' (2-7) 
^ i=i ' 

TT 2^x GM ( Xl 2xj -4-4 -2x1 + 3xi(x| + x§) nm -A 



+ 



R \ R 2R 2 2R? 

3G£ii 
~2R3 



where we only include the contribution from the quadrupole moment of each star, 
and neglect the higher multipoles. In this case, we may assume that each star of 
binary is an ellipsoid of its axial length a\ , 0,2 and . Index symbol Ai is that defined 
in the textbook of Chandrasekhar t3 and Aq = J2i^-i a i is calculated from 

poo fa 

a 2 a 2 a 3 I = = ajAo(a 2 , a 3 ), (2-9) 



where a 2 = a 2/ai and «3 = 03/^1 • Note that A{ is the function of a 2 and 03®. 
For ellipsoids, the pressure P is given by 



3 x 2 



^ = ^0 l-£3 • (2-10) 



\ 1 a i 1 

\ i=i 1 / 

Hereafter, we assume for simplicity that the internal motion of one star is the same 
as that of the other star, and we restrict its form given as 

u\ = — Ax 2 u 2 = Ax\. (2-11) 

a 2 a\ 

In this case, the vorticity u = 2Q + J2j,k e 3jkdjUk in the inertial frame is spatially 
constant everywhere inside the star as 

u = 2f2-h + ^)A^f2(2 + f R y, f R = - a l±A* (2 . 12) 
\ai a 2 J a x a 2 U 

Thus, the circulation in the equatorial plane becomes 

C = ira 1 a 2 fl{2 + f R ). (2-13) 

From the first tensor virial (TV) relation, we obtain the orbital angular velocity 

2GM lSGp_ 

From the second TV relation, we obtain the following three equations; 



A 2= I* _ 2 ^GpA l + ^ + Q 2 - 2^AQ, 
paf R 6 a\ 

pa^ R A a 2 



A 2 = — I - 2irGpA 2 - — r- + Q 2 - 2— An, 



2P „ GM 

0= ^- 2 ' GpA3 -->e- (245) 

Note that relations among R/ai, a 2 and 03 are derived from these three equations. 
The energy and the angular momentum are calculated as 



E= — 

5 



fi — —A \ a\+ (n-^A \ a 2 

«1 / V «2 J 



-^GpMA -^(M + 3 -±±y (2-16) 

j™tf + 4> a -^^ + M*n. (2-17) 

5 5 2 

In Eq. (2-16), the terms in brackets [■ • •] denote the spin kinetic energy of each star, 
the next two terms denote the kinetic energy of the orbital motion and the binding 
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energy by the self-gravity, and the final two terms denote the binding energy between 
two stars. The sum of the kinetic energy of the orbital motion and the binding energy 
between two stars mainly concerns the stability of the orbit and it becomes 

GM 2 3M*u 

Since a\ > 02,03 and as a result in > 0, the second term is positive definite. 
Thus, the energy has a minimum critical value E CT i t at a small radius due to the 
second term. As mentioned above, this just indicates that the tidal effect acts as the 
destabilization. 

Before presenting results on equilibrium states, we mention the possible state 
of internal motion of star. If the viscosity of the fluid is so small that the dissipa- 
tion time scale of the circulation by the viscosity (t V is) is much longer than i gw , the 
circulation conserves throughout the whole evolution. On the other hand, if the vis- 
cosity is so large that i v ; s is shorter than igw, BNSs will settle down to the corotating 
configuration. Kochanek, and Bildsten and Cutler showed that in order to achieve 
corotating binary dissipating the circulation, the viscosity must satisfy© 

^ 2o(R\ 2 ( m \Vl0km\ 5 / / / \ 

" z 10 UJ [mi;) {— ) g / cm/s - (219) 

where ao = (3M/47T/9) 1//3 and Mq denotes the solar mass. Since this value is much 
larger than the microscopic viscosity, E2P we should consider sequences of equilibrium 
states of constant C as realistic BNSs. 

By the way, for large R — > 00, fl — > 0, so that C — > — vr(af + af)A This implies 
that — A at R — > 00 is regardecLas the intrinsic spin angular velocity of each star (we 
denote it as i? s following LRS IMP ) . Since \f2 s \ for real NS's will be less than ~ 0.6 Qq, 
where i? = (GM/ag) 1 / 2 ,©!! l£7j should be less than ~ 4(GMo ) 1/2 . Note that 
for pulsars we have ever known, l3* \Q s \ ^ 0.2i?o- 

In Figs. 3, we show the energy and angular momentum as functions of Q/ Qq for 
C ' I ' (GMao) 1 / 2 = —2 ~ 2. We also show them for corotating binary for comparison. 
We see that there exists a critical angular velocity where the energy and angular 
momentum are simultaneously minima. As LRS showed, B the minima in the en- 
ergy and angular momentum along the constant-C sequence indicate the onset of 
dynamical instability. Thus, the minimum corresponds to the ISCO. As mentioned 
above, this minimum occurs due to the tidal interaction between two stars at small 
separation. Each star deforms and the quadrupole dependent term in the effective 
interaction potential between two stars (i.e., the third term in Eq. ( |2-3| )) becomes 
sufficiently large. 

Figures 3 also show some important features for the ISCO. First, the effect of 
the spin of each star is not so important in determining the ISCO as long as it is 
not very large, i.e., \C\/{GMa ) 1/2 < 1. Second, even for corotating sequence where 
the circulation does not conserve, there exists a critical angular velocity where the 
energy and angular momentum are simultaneously minima, and it locates near the 
minima for constant-C sequence. Although the minima for corotating sequence do 
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not mean the onset of the dynamical instability, but of the secular ins tabilityi (so 
we call it the innermost stable corotating circular orbit(ISCCO)), its existence is 
useful for analysis of the ISCO. This is because the angular velocity at the minima 
for corotating sequence is only slightly smaller than that for constant-C (small C) 
sequences and as a result, we can approximately find the location of the ISCO for 
constant-C sequences by investigating corotating sequence instead. 

2.2. Compressible case 

Next, we argue the effect of EOS. Since the location of the ISCO for binary 
stars is determined by the degree of the tidal deformation, it depends on the EOS 
of NS. NS becomes more centrally concentrated as the EOS is softer. Thus, the 
tidal effect is less important for softer EOS. This means that the effect of EOS is 
quantitatively very important in determining the location of the ISCO. Actually, 
using an approximate analytical model, LRS have shown that if we assume the 
polytropic EOS as 

P = Kp r ; T = l + -, (2-20) 

n 

where K and n are the polytropic constant and the polytropic index, -Rrsco is smaller 
for larger n, and for sufficiently large n 1), the ISCO disappears. 

Although LRS presented approximate locations of the ISCO for various n, we 
do not know the precise ones. Furthermore, their approximation is good only for 
small n, so that quantitative dependence on th e E OS may not be understood well 
for realistic NS whose EOS has 0.5 <; n ^ 1.0 To get precise locations of the 
ISCO, we need to perform numerical calculation, but we have not yet established 
an accurate numerical method of obtaining equilibrium states including arbitrary 
circulation. As we mentioned above, however, corotating sequence may be regarded 
as an approximate sequence of constant circulation if the spin of each star is not so 
large, and we have accurate methods of obtaining corotating sequences. ciP Thus, we 
here see the effect of the EOS by investigating corotating sequences. 

Equilibrium states of corotating binary stars are obtained by consistently solving 
the Newtonian potential U and the integrated Euler equation as 

K{n + l) P l ' n = U + ^-f2 2 ; w 2 =(^ Xl + ^j +x 2 2 . (2-21) 

In Figs. 4, we show the energy and the angular momentum as functions of the angular 
velocity for n = 0.25, 0.5, 0.75 and 1 (r = 5, 3, 7/3 and 2). The energy, angular 
momentum, and angular velocity are shown in units of GM 2 /ao, (GM 3 ao) 1//2 , and 
J?o = (GM/ag) 1 / 2 , respectively, where ao denotes the radius of the spherical poly- 
tropic star. In each figure, the circle of the largest fl denotes the configuration at 
contact of two stars. It is found from figs. 4 that with increase of n, fl at the energy 
and angular momentum minima increases and the minima tend to disappear for given 
ao and M. The angular velocity at the minima, ^N:iscco> is — (0.252 — 0.256) Qq 
for n = 0.25, ~ (0.265 - 0.269)12o for n = 0.5, ~ (0.279 - 0.283)^ for n = 0.75 and 
~ (0.296 - 0.298)j? for n = 1.0 Note that for incompressible case, it is ~ 0.248^ - 

*' Note that in the analysis by LRS,!' they somewhat overestimate J?N:iscco for compressible 
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From these numerical results, the angular velocity at the minima for compressible 
star of n = — 1 is found to be approximately expressed as (see Fig. 5) 

^N:iscco - 0.248(1 + 0.2n a )n , (2-22) 

where a ~ 1.5. Thus, the frequency of GWs corresponding to ^njscco is calculated 
from /n:ISCCO = ^NiISCCo/tt, which is 



M \ 1/2 /l5km\ 3/2 



Wco- 590(1 + .2n-)^— j ^— J Hz. (2-23) 

The above results concerning the dependence of the energy and angular mo- 
mentum minima on the EOS are for corotating binaries, but we may expect that 
similar dependence holds for constant-C sequences if the spin of each star is not so 
large. This is because change of the EOS mainly affects the density profile of each 
star almost independent of the spin and circulation. Actually, from the results by 
an approximate study of LRS,tf it is found that the following quantity is almost 
independent of the EOS; 

^N:isco for C = sequence (2-24) 



^N:iscco fo r corotating sequence 



Thus, we are allowed to conclude that even for constant-C sequences the angular 
frequency at the ISCO for compressible star will increase approximately with n L5 as 
Eq. (|2~22|). 



§3. First post Newtonian study: Incompressible case 

In §2, we briefly mention that the ISCO exists even in the Newtonian case if 
the EOS of NS is sufficiently stiff. In this section, we see how the location of the 
ISCO obtained in the Newtonian theory is affected by the 1PN correction of general 
relativity. We again make use of sequences of corotating binary stars. In this section, 
as a first step, we investigate equilibrium states for incompressible binary stars for 
which fairly accurate approximate solutions are obtained in an analytic manner (see 
details in Ref. 12)). 

3.1. Formulation 

Non-axisymmetric equilibrium configurations of uniformly rotating incompress- 
ible fluid in the 1PN approximation are obtained by solving the integrated form 
of the Euler equation and the Poisson equations for gravitational potentials consiSj- 
tently. The integrated form of the Euler equation was derived by ChandrasekharE3' 
and it can be written as 0)00 

J p c 2 J \ p J p 

star. For example, i?N:iscco = 0.279J?o for n = 0.5 and ,!?N:iscco = 0.314J?o for n = 1 in their 
analysis. 

*' Note that we use the standard PN gauge. 
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4 

TO 4 „ 4 



+ ^ (2t^C7 - X n + /^J i ^ + ^tf 4 + const, (3-1 



where Xq, Xq, and /?„ are the gravitational potentials in the 1PN order which are 
obtained by solving the following Poisson equations: 



APi 



/ P^ 
+ 



3P 



and 



+- 



^ * 2 
Z\P 2 = -4vrG/)X2 , 
AX = 4irGp(e + 2U + 
AX Q = 8irGpm 2 , 



XX + —J P 2 ,2 +a&P 1 ,l " Ui + - W 2 (Pl,2 + P^ 



(3-2) 



• (3-3) 



Note that in the incompressible = and p =constant. Thus, the left hand-side 

of Eq. (3-1) is 

p_±(pV 

p 2cAp 



(3-4) 



In the 1PN approximation, we have conserved quantities such as the cop served 
mass M*, the energy E, and the angular momentum J which are written asEZPlHP 



E 



p*d x 



M 1 + 



1 (v 



r v 2 i 

1 /5 



J 



+ - 



n 2 



i 



+ 3J7 W 3 x, 



P 



(3-5) 



CP ^,(3-6) 



1 



P 



c.-<i 1 + — I v* + 6U + e + - ] J> + 



d 6 x. 



(3-7) 



iW # conserves throughout the whole evolution of the system even if there exists a 
dissipation process such as the emission of GWs. Thus, in the case when we consider 
a sequence of a constant M*, it may be regarded as an evolution sequence of the 
system. 

Our purpose is to calculate the 1PN correction of the angular velocity, the en- 
ergy and the angular momentum for corotating binary. To get the solution in the 
Newtonian order, we solve Eq. (2-15) with A = using the angular velocity 



Q 2 



2GM 18G^n 
+ 



P 3 



P 5 



(7 2 



(3-8) 
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In the 1PN approximation, we can expect that the following types of quantities will 
be the main terms of the angular velocity in the 1PN order: 

GM GM GM GM GMal GM , GMal GM , . 

where ag is a typical radius of the star and we use the relation J-n ~ Moq. We will 
derive the details of four types of terms shown above in the correction of the energy 
and the angular momentum neglecting the effect of higher multipole terms. 

As shown in Ref. 12), the neglection of the higher multipole terms brings us a 
great advantage so that we can regard each star of the binary as an ellipsoid whose 
shape is the same as that in the Newtonian order. This is because the deformation 
is induced by higher order effect of R~ . Hence, we hereafter perform calculation 
setting density profile of each star as an ellipsoid of its axial length a%, a 2 and 03. 

3.2. Gravitational potentials 

As we did for the Newtonian potential U, we decompose PN gravitational po- 
tentials into two parts; one is the contribution from star 1 and the other is from 
star 2. In the following, we denote the former part as (j) 1 , and the latter one as 
(j) 2 ^ 1 , where 4> denotes one of the potentials. We also define cf) 1 ^ 2 and 2 ^ 2 as the 
contribution from star 1 to 2 and star 2 itself, respectively. 

3.2.1. X 

As we mentioned above, the 1PN potentials are divided into two parts as Xq = 
Xq^ 1 + Xq^ 1 , and we consider them separately. 
• Contribution from star 1 

Aq^ 1 is derived from the Poisson equation 0) 



AX^ 1 = AirGp 



1 




2vGp K-E |+-(l-E^I^^ 



(3-10) 



and the solution becomes 



X^ 1 = -aoU 1 - 1 + a x D 1 + £ m D u - — ( 2D n - D 22 - D 33 



M 
R± 



2D ni + 3£> 122 + 3D 133 ), (3-11) 



where 



„ , 3P 2GM 3G+n 
a = 2ttG P A + — ^ + — — + — 3-12 
p R R a 

2GM 9Gfn 
«l = ^^ + — S3—. ( 3 ' 13 ) 



R 2 i? 4 ' 



op 

m = 2irGpAi + -4. (3-14) 



12 



Di, Da, and Dm are the solutions of equations 

ADi = -AnGpxi, 
AD U = -AnGpxl 
ADiu = -inGpxix 2 , 

and the solutions at star 1 are© 

A = nGpaj [A,, - £ A u xf 



(3-15) 
(3-16) 
(3-17) 

(3-18) 



Da = irGp 



An — Amxf j x 2 



+-a 2 ( Bi - 2 ^ A/2; 2 + E E B Umxfx 2 m 



l m 



(3-19) 



An = nGp 



5 Ain-^im^ ^ 



- 2 E + E E Bnimxfxl) x x , (3-20) 



Dm = irGp 



\a\ ( ^lii - E Maixj 



X\X; 



+ 7 a l a i ( B li ~ 2 E ^liJ^i + E E B lUrnxfx 2 m 



I m 



where Ajfc- - are index symbols defined by Chandrasekhar 
Contribution from star 2 
The equation for Xq^ 1 is 



(3-21) 



^■T 1 = ^Gp 2nGp (ao - E A W ^j + ^ (l - E | ) + 2C/1 ^ 



(3-22) 



where yi = — (xi + R), y2 = X2 and 2/3 = X3. The solution is expressed as 



1 



where D 2 ~". 1 are calculated from the same equations as the case of Dij...,c3 i.e., 



-AirGpytyj 



M I 9 n2-.l n 2^1 n 2^1 
^1^11 _LI 22 -^33 )' 

P i.e 
(3-24) 



The solutions of Dfr* are 



2—1 



i? 2 



1 1 / , 2xi + 6xf - 3aj - 3aj + ^- 3 -, 



2i? 2 



(3-25) 
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n2-»i Glut x\ 2x\ — x\ — x\ — 2x\ + 3xi(x 2 + x|) ^/n-^l 
D « = ~R~ 2P^ + 2R? + 0{R ] 



+ 2P 3 



where 



(l-^ + 0(W), (3-27) 
^-iiii = ^iiii ~ o X! I™ 11 ' (3-28) 



3.2.2. X fi 

• Contribution from star 1 



The equation for Xjf* 1 is 



AXjp 1 = %irGp(x\ + xl + Rx 1 + ^ r ). (3-29) 
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Then the solution is 



XjT 1 = -2 ^ n + D 22 + + ^-U 1 - 1 j . (3-30) 

• Contribution from star 2 

The equation Xf^ 1 may be written as 

AXfT 1 = 8^G P ^ + y| - P yi + ^ . (3-31) 

Then, using Dfj^ 1 , the solution is easily derived as 

xt 1 = -2 (^r 1 + ^V 1 - rdi- 1 + ^V 2 - 1 ) , (3-32) 

R 2 tt2 ^ 2GI n 2GI 22 (, x x 



= - T ^-i^-^ I--)- (3-33) 



3.2.3. /3 V 

$ v is obtained from Pi and P 2 which satisfy, 



AP t = -AttGp ( xi + , (3-34) 



2 

Z\P 2 = -AitGpx 2 . (3-35) 



Pi is also written as P/ - * 1 + Pj 2_>1 , where 



P^^-Di + I^ 1 , (3-36) 
P^ 1 = D 2 (3-37) 
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and 



1 1 2 



■2->l 



2-»l 



(3-38) 
(3-39) 



3.3. The Post- Newtonian angular velocity and definition of mass and center of mass 

As in the Newtonian case, the orbital angular velocity in the 1PN order is also 
derived using the first TV equation as 







— d 3 x = / pU^x + -MQ 2 + - pX !d 3 x 
ox\ J 2 & J 

J p^2w 2 U tl +4 Xl U + 2RU - X n ,i + 9 ,^ja9x 



a R, s MR 3 

+^ -(3J 11 +/ 22 ) + — — 



(3-40) 



where ", fc" denotes the partial derivative with respect to xt- Since all the gravita- 
tional potentials can be written as the polynomial form of x^, we can immediately 
perform the integral shown above and obtain Q 2 as 



Q 2 = 



2GM 



1 + | UnpAo - ^ - ^3 (28a? - 14a* - 9a 2 ) + 0(RT 



R 5 



c 2 [ 
1 + 



28 

15? 



TrGpA + 0(R- 2 )^j , 



(3-41) 



where we make use of the relations A = J2i-^-l a f an d Eqs. (2-15) with A = to 
simplify the expression. 

Next, we discuss the mass of the system and definition of the center of mass 
for each star. This is because there are several definitions of them in the 1PN 
approximation, and we should clarify the difference between the similar ones. 

First, we consider the conserved mass which is defined as 



= 1* 



xp 



= M 



1 + 



1 + 7>{Y + 3U ) 

G ( 13M 127rpA 



^ AR 



+ 



+ 2W (34a '" lla 2" 15a ' ) + ° ( ^ 



(3-42) 



Using M*, Q 2 becomes 



Q 1 



2GM* 



n G 

c z 



R 3 
18G(*n), 



2irpA 11M* 2 2 2 1 ^, i? _4a 
5 ^ - ^3(90ai -39a 2 -33a 3 )| + 0(i? ) 



i? 5 



13M* 
4i? 



+ o(ir 2 ) 



(3-43) 
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where (J-n)* = M m -t\\/M. Thus, Q 2 looks as if it depends on the internal structure 
of the star even in the limit of ai/R — > 0. Since we believe that in the EOM for the 
point particle, the quantities depending on the internal structure does not appear, 



M* is not desirable to describe the EOM for the point partic 
when the EOM is derived, one usually adopts the PPN mass 



Instead, in the case 
2 3 defined as 



MppN 



d 3 xp 



1 v 2 



1. 



M 



1 + 



1 + ^[Y + 3U ~2 Useli + 

G ( 13M _ . M 



y self 



\ 4R 



+ 2tt P A + ^ (38a 2 - 7a 2 - 15a 2 ) + 0(iT 5 ) 

(3-44) 



where U se \{ and v se \f are the self-gravity part of the Newtonian potential and the spin 
velocity of each star. When we rewrite Eq. (3-41) using the PPN mass, the orbital 
angular velocity does not depend on the internal structure of the star and it agrees 
with that for the point particle® in the limit of ai/R — > as 



9 2 



2GAfppN 



R? 



1 + 



G 



llMi 



PPN 



Mi 



PPN 



2R 



20R 3 



(94af - 35a^ - 33a^) + OfTT 4 ) 



+ 



18G( j : n)ppN 
R 5 



1 + 



G 



2 13M PPN (0 / E) _2> 

■^ P A —^- + O(R , 



(3-45) 



where (J-h)ppn = AfpPN^n/-^'- Thus, when we compare the present results with 
the point particle calculations, we should use the PPN mass. In the present case, 
however, Mppn is not a conserved quantity although M* is. When we consider a 
sequence of equilibrium configurations as an evolutionary sequence, we should fix 
M*. 

Next, we consider the definition of the center of mass for each star. In the PPN 
formalism, it is defined ascS 



x PPN 



1 



MppN 



d xpx 1 



1 



1. 



1 + ^7 -rr + 3C/--i7 self + 



y self 



(3-46) 



x\ coordinate of the center of mass for star 1 deviates from to 



G ( 2Ma\ 
c 2 [~~5R 2 ~ 



+ 0(R~ 



(3-47) 



Thus, in the PPN formalism, the following orbital separation should be used; 



Rppn = R 



1 + 



G f AMa 2 
1? I 



+ Q(R 



-5^ 



(3-48) 



It is worth noting that when we define the center of mass by the conserved mass 



as 



JT I ,f ' xp 



(3-49) 
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the result is the same up to 0(iT 4 ). Thus, in this paper, we do not have to distin- 
guish from -Rppn- Even in the general case, the difference between R* and -Rppn 
is expected to be small. 

Using R* and/or -Rppn, ^ 2 is rewritten as 



o 2 



2GM* 
Rl 



1 + 



G J 2tt P A 11M» 



2R* 



2QRI 



138a 2 - 39a 2 . - 33o§) + 0(R; 



+ 



18G(-t 



iu* 



G 



13M* 
4i?* 



(3-50) 



or 



J? 



2 _ 2GMppn 
n ppN 



1 + 



2i? PPN 

_ -^PPN 
20i?pp^ 



;i42a 2 - 35a 2 . - 33a§) + 0(U£| N ) 



18G(* 



HjPPN 



-Rp PN 



, g r 2 , 



13Mi 



PPN 



4-Rppn 



0(i2 



-2 i 
PPN; 



(3-51) 



Here, we should note that the effect of the spin-orbit coupling terms will appear in 
Q 2 from 0(i?pp N ).EP According to the 1PN study, it becomes 



Q 2 



2GAfppN 

-rt ppN 



1 



2GiWppN 

r 2 E>3 



(a? + ^ 



(3-52) 



where we omit other terms which are not related to this discussion. Eq. ( |3-52| ) 
shows that the terms of 0(i?pp N ) in Eq. (3-51) cannot be explained only by the 
spin-orbit coupling term. This will mean that there appears a new effect, say the 
1PN quadrupole one, in Eq. (3-51). El 

3.4. The energy and the angular momentum 
3.4.1. The Total Energy 

The 1PN total energy is calculated from E = E& + Ep^/c 2 , where 



zpdef 



TP' 



def 
PN 



P 



P 



1. 



U d 3 x, 



1 

—v 
2 

V + 5 2 P 2 _ 5 2 l- 2 

8 2 p 2 2^ 



(3-53) 
(3-54) 



*' In this paper, we use the PPN mass in order to compare the angular velocity with that of 
the point particle case. However, the PPN mass is useful only for the point particle case because it 
does not conserve if we take into account the tidal forces on each star of binary. In the context of 
this paper, we have to use the mass which is valid even if the tidal forces exist. Unfortunately, no 
one has proposed such a mass as far as we know. 
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For the 1PN corotating binary, they become 

4ttG P A Q 2 , 2 2x r2 ^ 2 GM 3G4- 



E N ei = M 



= 2G 2 M 



5 +T (ai + fl2)+ 4 
34, . . 2 llMnpAo 7M 2 



ii 



R 



R? 



+ 0{R~ 



+ 



+ 



3R 



MTTpAp J 68 2 2 N 
M 2 



32R 2 
61 

7 M 



240P 4 



(302af + 59a2 - 209a|) + O(ir ) 



(3-55) 



If we substitute Q 2 into the above formulas, E may be rewritten as En + -Epn/c 2 
where 



-En = M 



^GpA, GM 3G1 U f2 2 2 2 
5 -2-R + ^R^ + ~ {ai + a2) + 0{R 



E PN = 2G 2 M 



34, 2 19MnpA 25M 2 



21 



6R 



32R 2 



, MttpAq (22 2 . 2^ 158^11 
+ -R^l21 (ai + a2) -^17 



+ 



M 2 



240i? 4 



(26a? + 35a 2 . - 155a 2 ;) + Q(R 



-5\ 



,(3-56) 



(3-57) 



3.4.2. The Total Angular Momentum 

We can calculate the 1PN total angular momentum from J — <^n 6 ^ ^pnA^' 



where BE 



T dcf 



T dof 



P 



v v [ v 2 + 6U + - ) + ^ 



P 



erx, 



(3-58) 
(3-59) 



and v v = Qvo 2 . For the PN corotating binary, they become 



T dcf 



164 



R\pA + 5i?Af + — TTpA (al + a 2 ,) 
1U5 

+y^(20a 2 + 15a 2 - 11a 2 ) + 0(R- 



(3-60) 
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If we substitute J? 2 into the above formulas, J may be rewritten as Jn + Jpn/c 2 
where 



, r ~ ( R 2 aft + at 
J N = 2Mf2 N — + 1 2 



(3-61) 



Jpn = GMQ~n 



\r\ P A q + -^-RM + (2018a 2 + 2081a 2 . + 21a§) 
2 lb 1050 

+-^-(122a 2 + 85a 2 2 - 97a!) + 0(iT 2 )l . (3-62) 
80x1 



3.5. Construction of equilibrium sequences 

To construct equilibrium sequences of 1PN corotating binary stars fixing M* 
and p, we use the following method: 

Step (1): Using Eqs. (2-15) with A = 0, we numerically calculate equilibrium se- 
quences in the Newtonian order. Up to this stage, a 2 , 0:3, and R = R/a± are 
determined. 

Step (2): a\ is determined from the condition M* =constant using Eq. (3-42): 



3M* 



4'7rpa2Q ? 3 



■KGp ( 3M* 



2/3 



I 4npa 2 a 3 J 1 



/ 12A 13a 2 a3 . «2«3 ^, ni 2 K 

= ^-(34-lla 2 -15a 3 ) 

5 3R 15R 3 J 

(3-63)' 



Step (3): After substituting Eq. (KMS3D into Eqs. (3-41), (3-55), and (3-60), we 
rewrite the 1PN expressions for the orbital angular velocity, the energy and the 
angular momentum as 



Q 2 
E 
J 



- o | GM* h2 



o 2 

2 J^PN) 



era 



E 



(GM2/a») 
J 

(GM*a*)V 2 



-En H — 5 — Epn, 



c 2 a* 



Jn H — 5— - Jpn, 



(3-64) 
(3-65) 
(3-66) 



where 



a* 



ni 



o 2 



/3M,y/ 3 



V 47rp 



J 



= a 203 



2 6 ,„ 9 

H ^(2 - a - a 

R 3 5R^ 2 



GAL 



at 

(a 2 a 3 ) 4/3 



3 A 







9 , 42 A 2 2 

+ tt^t (2 - a 2 - "3 



R 3 a 2 a 3 2R 4 25R 5 a 2 a 3 
1 



5i? 6 



,-(28 - 14a| - 9ag 



(3-67) 
(3-68) 
(3-69) 



(3-70) 
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E N = (a 2 a 3 ) 1/3 



3 A, 



5 a 2 a 3 2R WR 3 



+^Ml + a|) 



£pN = («2a3) 2/3 



5«2«3 

3 / jp ' 
140 \ 0:203 



55 1 
+ — -77 + 



48i? 2 700.R 3 "203 



^(194 + 215a^ + 165a?) 



J N = 2(0203) 



-1/6 



O 2 
o 2 o 3 



.1/2 




(3-71) 
(398 + 401o^ + a§) 
(3-72) 

(3-73) 



J PN = (o 2 03) 1/6 
+ 1 4) 



o 2 


1/2 - 


0203 





3~n A 83 - 

— R 2 — — + — R 

8 0203 48 



1400 0203 



(338 + 401a^ + 21a$) 



240i? 



(494 + 155a 2 . + 141a§) 



(3-74) 



Then, using the numerical value of 02, 03, and R determined at step (1), we calculate 
the sequence of the angular velocity, the energy and the angular momentum as 
functions of the orbital separation. 

We repeat this procedure changing the mean radius of each star a*. Once a 
sequence is obtained, we search the minimum point of the energy. If we find it, we 
call it the ISCCO. The 1PN approximation is valid only for small C s = GM*/c 2 a*, 
i.e., the characteristic value of the compactness of each star. Thus, we hereafter show 
results only for small C s . 

In Figs.6(a) and (b), we show E = E/(GM 2 /a*) and J = J / (GM^a*) 1 / 2 as 
functions of Q/Q*. The figures show the important fact that the angular velocity at 
the ISCCO increases approximately in proportion to C s . 

In Fig.7, we show Q = Q/Q* at the ISCCO as a function of C s . We show Q at 
the energy minimum as well as that at the angular momentum minimum. The figure 
indicates that two minima are almost coincident, but slightly different. We guess 
that this disagreement is due to the neglect of the effect from the higher multipole 
deformation. In any case, we may expect that the ISCCO locates near two minima. 
Fig.7 clearly shows that the orbital angular velocity at the ISCCO increases almost 
linearly with increase of C s as 



^ISCCO = ^N:ISCCO (l + C PN (n)a 



(3-75) 



where Cp-^(n) is a constant which depends on the EOS (see next section) and for 
the incompressible case, Cpn(0) ~ 0.5. 

We briefly argue the meaning of the result for i^iscco using Eq. (3-72). In 
Eq. (3-72), there appear mainly two kinds of effects in the first and second terms, 
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respectively: The first term concerns the internal structure of each star and denotes 
that each star is forced to be compact due to the 1PN gravity because of its negative 
definite character. Thus, the first term acts as the stabilization because the tidal 
effect is less important in the case when the star becomes compact. On the other 
hand, the second term concerns the gravity acting between two stars because it exists 
in the limit ai, 02, 03 — > 0. This term tends to act as destabilization of circular orbits 
because of its dependence of R (0(R~ 2 )) and its positive definite character. In this 
way, there exist two opposite effects in the 1PN corrections, and the present result 
shows that the former effect dominates over the latter one. 



§4. First post Newtonian study: Compressible case 



For compressible stars, we can no longer solve the integrated form of the Eu- 
ler equation analytically although we make approximate solutions assuming density 
profile of each star appropriately. H2P e3 To obtain accurate equilibrium states, we 
need numerical calculation. In this section, we show numerical results for sequences 
of 1PN corotating equilibrium binaries obtained by solving the integrated form of 
the Euler equation as well as the Poisson equations for the gravitational potentials. 
We use a numerical method which is developed in previous papers. IHf 

In the compressible case, the integrated form of the Euler equation for uniformly 
rotating fluid is 



K ( n+1)p Vn _ J_(x in+l)p l/ny 



= C/ -fr+j^-+^ ( 2w2u ~ X + &>) + |J^ 4 + const, (4-1) 

where we assume the polytropic EOS, i.e., 

P = Kp W/n and £= r ^ ( 4 . 2 ) 

P 

Hereafter, we set n = 0.5, 0.75 and 1 because real NSs will have n = 0.5 ~ l-H* 
The gravitational potential in the 1PN order are obtained by solving Eqs. (3-2). 
The conserved mass, the energy, and the angular momentum are defined as Eqs. 
(3-5)-(3-7). Following a previous section, we hereafter show the angular velocity, the 
energy and the angular momentum in units of I?* = (GM t /o]) 1/ ' 2 , GM^/a* and 
(GMfa*) 1//2 , where 

for n = 0.5, 



for n = 0.75, 




0.9960G J 

1/2 




for n = 1, 
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(4-3) 



We also define the center of mass for each star as 



Xl * = W I p * Xld3x > ( 4 ' 4 ) 

and as a result, the orbital separation is defined to be R* = 2x\*. In numerical 
computation, sequences of equilibrium states are obtained fixing M* and K. 

First, we demonstrate that our numerical results are reliable by showing the 
relation between the orbital separation and the angular velocity. If we neglect the 
contribution from the quadrupole deformation of each star, the angular velocity in 
the 1PN approximation can be approximately expressed asEJ E!P 
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IGM* 

m 



l G ( 3 - n M* 11M* 



c 2 I 2(5 — n) a* 8R* 



(4-5) 



Numerical results should agree well with this analytic formula. In Fig. 8, we show 
the relation between R*/a* and 17/ '17* for n = 0.5 and C s = GM*/a*c 2 = (open 
circles), 0.0327 (filled circles) and 0.0654 (open squares) as examples. Dotted lines 
are drawn using Eq. ( f4-5| ) as counter parts of the numerical results. The figure 
shows that numerical results agree with the analytic formula well except for slight 
deviation caused by the deformation of each star. 

In Figs. 9, we show the energy and angular momentum as functions of 17/17* for 
sequences of n = 0.5(a), 0.75(b) and 1(c). In Fig. 9(a), open circles, filled circles, 
and open squares denotes sequences for C s = 0, 0.0327 and 0.0654, in Fig. 9(b), they 
denote sequences for C s = 0, 0.02 and 0.04, and in Fig. 9(c), they denote sequences 
for C s = 0, 0.0167 and 0.0333, respectively. It is shown that the minima of the 
energy and the angular momentum almost coincide, and 17/17,,, at those minima 
(i.e., ISCCO) increases with increase of C s as in the incompressible case. 

In figs. 10, we show the orbital angular velocity at the ISCCO as a function of 
C s for n = 0.5(a), 0.75(b) and 1(c). It is found that as in the case n = 0, l7iscco 
almost linearly increases with increase of C s as 



^ISCCO = ^N:ISCCO (l + C PN (n)C s ), (4-6) 



where Cpn(0.5) ~ 1.1, Cpn(0.75) ~ 1.8, and Cpn(1) ~ 2.5. Thus, not only 17n : isccO) 
but also Cpn changes with the EOS, and for softer EOS, 17iscco sensitively increases 
with C s because the 1PN effect of the self gravity is larger for the softer EOS. 
The summary of §3 and §4 is as follows. 

1. 17iscco increases with the increase of the 1PN GR correction. 

This is because each star of the binary becomes compact due to the strong GR 
self-gravity and as a result, the tidal effect becomes less effective. 

2. The feature mentioned at 1 is more remarkable for larger n, i.e., for softer EOS. 
The reason for this behavior is that for softer EOS, each star has the larger 
central density and is susceptible to the GR correction. 
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These conclusions are obtained by analyzing equilibrium sequences of corotating 
binary. However, it seems that these features may not be very sensitive to the spin 
and circulation of each star because the above properties mainly originate from the 
change of the density profile due to the GR correction. Therefore, we may expect 
that the similar trend will hold for constant-C sequences although more investigation 
is needed for this point. 

§5. GR orbital effects: Pseudo-Newtonian potential approach 

In previous two sections, we investigated the stability of corotating binary in 
the 1PN approximation, and showed ^iscco increases with the increase of the GR 
correction. This may be viewed from another point. In the 1PN approximation 
ISCO does not exist in the point particle limit. E^P Therefore with the increase of 
relativity parameter C s , i?isco should increase since the increase of the relativity 
parameter means the way toward the point particle limit where the circular orbit 
exists for any radius. On the other hand, if we include the GR effect to the gravity 
acting between two stars in which ISCO exists even in the point particle limit but we 
do not include the self-gravity of each star, i7isco decreases compared with .i?isco 
in the point particle limit. This seems to apparently disagree with the results Jn 
previous sections. In the following, we first review the result of the work by TN,lliP 
which is one of such studies. (We note that there is the other work by LW,© in 
which the qualitative features of the results are essentially the same as those by TN.) 
In the last subsection, we clarify the reason for this apparent disagreement. We will 
show that this disagreement simply comes from the difference of two treatments. We 
also suggest that if 1PN effect of the self-gravity is included in the formalism of TN, 
■^isco w iH increase from the results without the 1PN self- gravity effect. 

5.1. Outline of method 

Here, we consider the so called Roche-Riemann problem© as a model of a 
binary that consists of a finite size star (star 1) and a point-like gravity source (star 
2). To mimic the GR potential between two stars in which circular orbits are unstable 
for small radii, we use the pseudo-Newtonian potential proposed by Paczyrisky and 
Wiita.E§ Details of the method by TN are as follows. 

First, we regard star 2 as a point-like star of mass mi and denote the gravi- 
tational potential by it as V^r). To include the effects of general relativity phe- 
nomenologically, we adopt pseudo-Newtonian potential as Viir) (see below). Then, 
as star 1, we adopt the incompressible and homogeneous ellipsoid of its axial length 
ai, a>2 and 03, mass m\ and the density p±. The structure of star 1 is determined by 
the Newtonian hydrostatic equation. Thus, the self-gravity part of the gravitational 
potential V% is the same as that in Eq. (2-7). 

5.2. Basic equations 

We use the tensor virial method to derive the equations necessary for deter- 
mining the orbital angular velocity and for constructing equilibrium configurations 
of star 1. We choose the coordinate system such that the origin is at the center of 
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mass of star 1 and the xi-axis points to the center of mass of star 2 which is located 
at (— R, 0, 0). The X3-axis is chosen to be the axis of the orbital rotation. In the 
frame of reference rotating with i?, the hydrostatic equation of star 1 is written as 



dm 



dP 







dxj Pl dxi 



Vi + V 2 + 



Q 2 



rri2R 
m i + m 2 



+ Xl\ + x 2 



i=i 



(5-1) 



where m, P, and R are the internal velocity, the pressure, and the separation between 
two stars, respectively. 

Following Chandrasekhar, t3 we expand the interaction potential V 2 in the power 
series of x^ up to the second order assuming R 3> a%,a 2 and a 3 . We assume that 
the potential V 2 depends only on the distance r from the center of mass of star 2 as 



where r is given by 



V 2 = V 2 (r), 



[{R + Xl f + x 2 2 + xl) 1/2 . 



The expansion of V 2 (r) near the coordinate origin becomes 



V 2 = (V 2 ) + 



dV 2 
dr 



x\ + 



1 d 2 V 2 
2\ dr 2 



xi + 



1 (dV 2 
2R 1 dr 



{x\ + x 2 ), 



(5-2) 
(5-3) 

(5-4) 



where the subscript denotes the derivatives at the origin of the coordinate. From 
the force balance at the center, we obtain the orbital angular velocity as 



m 2 R 



-Q 2 



dV 2 
dr 



(5-5) 



m\+m 2 

where 5 is the quadrupole term of the interaction potential!) . If we take the New- 
tonian potential 

Gm 2 



V 2 (r) = 

as an interaction potential, 5 is written as 



(5-6) 



WR 2 



(2a? - al - a\ 



(5-7) 



Substituting Eqs. (5-4) and (5-5) into Eq. (5T), we have 



dui 



dP 







j=l 3 



+ 



dxi dxj 

1 (dV 2 
2R\ dr 



Vi-S 



if 

\ or 



2 2^ ^(d 2 V 2 

Xl + -n 2 ( x 2 + x 2 ) + -\ — 



r 2 
1 



x 2 + ^i) 



+ 2pif2^2eu 3 ui. 
i=i 



/ o 
(5-8) 
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Multiplying xj to Eq. (5-8) and integrating over the volume of star 1, we have 
= 2T ij + W ij + L 2 +(^\ U^. + l^ + I^ \s 2i I 2j 



1 ( \ ^ f 



x + s i:i n, (5-9) 



where 



Tij = — J piUiUjd 3 x : Kinetic energy tensor, (5 - 10) 

Wij = / pi^-^-Xjd 3 x : Potential energy tensor, (5'H) 

-U*. — ,,, 

n= [ Pd 3 x. (5-13) 



and 



In Eq. (5-9) there appear no terms related to 5. Since it is possible to take the 
coordinate system moving with a constant velocity, the term proportional to 5 in Eq. 
(5-8) can be vanished. Equation (5-9) is the basic equation to construct equilibrium 
configurations of star 1. Note that in this model, the GR effects generated by star 2 
are included in the orbital angular velocity and the tidal potential at star 1. 

Using Eqs. (2-7), (|2-10| ), and ( 2-11 ) for the self-gravity part of the potential V±, 



the pressure P, and the internal velocity Uj, respectively, we obtain three equations 
from Eq. (5-9); 

-A" = ^ - 2irG Pl A l +0 2 - 2^ An, 

ptaf \ dr Jo fl i 

-A' = ^ - 2nG Pl A 2 + + " 2 -^< 

= ^-2^^3 + 1^] . (5-14) 
Pi a 3 R \ dr L 



If we take the Newtonian potential as an interaction potential, Eq. (5-14) reduces 
to Eq. (2-15). Combining these equations with the angular velocity equation 

where p = mi/m2, we can construct equilibrium sequences. 
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5.3. Modified pseudo-Newtonian potential 

Although there are a wide variety of choices of V 2 {r) to incorporate GR effects 
phenomenologically, we use the so called pseudo-Newtonian potential modifying the 
original form proposed by Paczyhsky and Wiita. © This potential fits the effective 
potential for test particles orbiting Schwarzschild black hole quite well, as we will 
show later. We will use the modified pseudo-Newtonian potential defined by 

V 2 (r) = _ Gm2 , (5-16) 

T '"pseudo 

rpseudo = r s {l + g(p)}, (5-17) 
_ 7A9p 10 V , 29.3p 3 

m ~ 6(l+p)2 3(l+p)4 + 6(l+p)6' 1 j 

_ 2GM tot 

r S = 5 ' ( 5 ' 19 ) 

M tot = m 1 + m 2 , (5-20) 

where g(p) is a correction term to fit the ISCOs obtained from the hybrid 2PN EOM 
of Kidder, Will, and Wiseman c3 . For p = 0, the modified pseudo-Newtonian po- 
tential agrees with the original pseudo-Newtonian potential proposed by Paczyhsky 
and Wiita. 

Figure 11(a) shows effective potentials (solid lines) and locations of circular orbits 
(dots) in the modified pseudo-Newtonian potential (p = and r pseu( j = r s ) and in 
the Schwarzschild metric. Although by this choice of the parameter (r pseu do = r s), 
the locations of the ISCOs in the modified pseudo-Newtonian potential agree with 
those in the Schwarzschild metric, the angular momenta at the ISCO are different, 
that is, the angular momentum in the modified pseudo-Newtonian potential («/ pse udo) 
for p = is (9/8) 1 / 2 times larger than that in the Schwarzschild metric (Jsch) at 
the ISCO. Therefore in Figs. 11(a) and (b) we compare circular orbits with different 
angular momenta related as 




Rpseudo — o <^Sch- 0>"21) 



From Fig. 11(b) we see that the radii of circular orbits of the modified pseudo- 
Newtonian potential agree with those of the effective potential around Schwarzschild 
black hole within 10% accuracy near the ISCO. This is the reason we believe that 
the modified pseudo-Newtonian potential defined here expresses the effect of general 
relativity within 10% accuracyor so. 

In the pseudo-Newtonian case, the quadrupole term 5 is written as 



10 



2 (3R Rpseudo )\R r pseudo) r 2 , 2\ 
2o a ^ (°2 + a 3/ 



\R Rpseudo) 



1 

.(5-22) 



When we take the limit r pseu d — > 0, Eq. (5-22) recovers the Newtonian potential 
case (Eq. (5-7)). 
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5.4. The energy and the angular momentum 

In this subsection, we show the total energy and angular momentum. We regard 
the minimum point of the total energy or the total angular momentum as the location 
of the ISCO. 

The total energy is formally written as 



E = T + W + Wi, 



(5-23) 



where T, W and Wi denote the kinetic, the self-gravity and the interaction en- 
ergy, respectively. In the pseudo-Newtonian potential approach, T and W have the 
same form as that in the "pure" Newtonian case, but the interaction energy has the 
different form as 



Wi = - 



Gmim,2 

1' Rpseudo 

2(R- r pseudo ) 3 



2/ n - 



R Rpseudo 

R 



(I22 + -J33) 



Thus, the total energy is 



10 



ai 



a l + (f2-^A 



02 



+ 



miR 2 



Gm\m2 



Gmo 



-ft Rpseudo 

where the angular velocity is written as 

G(mi + m 2 ) 



2/ n 



2(1 +p) 

R Rpseudo 

R 



(5-24) 



-irGpimiAo 
5 



(I22 + hs) 



(5-25) 



1 + 



W(R - r pseudo ) 



R{R Rpseudo) 

(3i? — rp Seudo ) (i? — rp Seudo ) 2 . 2\ 
3^ ( ° 2 + a?>) 



(5-26) 



The form of the total angular momentum in the pseudo-Newtonian potential 
approach is the same as that in the pure Newtonian case as 



m\R 2 mi, 2 2 . 2 

J = ii H — —{a-, + an)S2 — -mia^A. 

l+p 5 5 



(5-27) 



Note that the pseudo-Newtonian effect is included implicitly in the angular velocity. 
5.5. Results 

In previous subsections, we described equations needed to get equilibrium se- 
quences and to determine the location of the ISCO. Here, we show the results for 
the equal mass binary (p = 1) in the corotating case (/# = 0) and the irrotational 
case (f R = -2). 

In Figs. 12 and 13, we show E(&) and J(b) as functions of I? in the corotating 
and irrotational cases, respectively. Note that the minima for the energy and angular 
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momentum denote the ISCCO for f R = case, and the ISCO for = —2 case. In 
Fig.14, we show fl = f2(Gmi/c 3 ) at the ISCO and/or ISCCO as a function of the 
compactness parameter of star 1 Gmi/aoc 2 for p = 1. Note that this normalization 
is different from that in Fig. 7 and so on. In this normalization, the orbital angular 
velocity of the Newtonian order (the value at Gmi/aoc 2 = 0) becomes zero, and 
it is not convenient to compare Q with the Newtonian angular velocity. However, 
the approach of TN is based on the point particle binary case, i.e., the case of 
Gmi/aQC 2 — > oo. Then, it is natural to choose the normalization i? in this section. 

From these figures, we find that due to the finite size effect of the star, Q at the 
ISCO (ISCCO) is always smaller than the point particle limit (17=0.0402) i.e. the 
finite size effect of the star always destabilize the binary. However as a function of 
the relativity parameter, fl at the ISCO (ISCCO) increases which is qualitatively 
the same result as in 1PN treatment. In this sense there is no disagreement between 
1PN and the result in this section apparently. 

5.6. 1PN calculations without 1PN self-gravity effect 

One of the unsatisfactory points in TN is that the PN effect of the self-gravity is 
not included. While in 1PN calculations in §3 and §4, 1PN effect of the self-gravity 
as well as the gravity between two stars are consistently taken into account. In 
§4 we suggested that 1PN effect of the self-gravity will have stabilizing effect. To 
expect PN effect of the self-gravity in TN in this subsection we will calculate i? at 
the ISCCO removing the 1PN self-gravity terms in the calculations of §3. 

5.6.1. Equations 

As mentioned above, we only incorporate the PN effect only to the gravity acting 
between two stars in this approximation. So that equations we use are quite similar 
to those in §2.1 except for the angular velocity which we use 



~ 2GM 



l + ^2(2a?-«2-al) 



9G 2 M 2 

W (5 ' 28) 



Here the PN effect appears in the last term. (Note that for a±, a-i, 03 — > 0, the angular 
velocity agrees with that for two point masses.) As in §2.1, we use Eqs. (2-15) to 
determine «2 = 02/01 and 03 = a^/a\. 

In calculating equilibrium sequences, we must fix the conserved mass M* (see 
Eq. (3-42)). In this section, the conserved mass should be written as 

*-"( 1 + w> 

where we neglect the stellar structure terms in the 1PN order. 

As was done in previous sections, we regard the minimum of the energy as the 
location of the ISCO. In this approximation, the energy is written as 

3GM 2 „ GM 2 GM 2 . 2 2 2 , 
E =-^aT A0 -^R + T0W^- a2 -^ 
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M 
~5~ 



{al + a 2 2 ){ 1 + 




4afa2 



25G 2 M 3 
16i?V : 



(5-30) 



where ao denotes the radius of the spherical star in the Newtonian order. 
5.6.2. Results 

In the following calculation, we normalize the angular velocity and the energy 
by (GM,/^) 1 / 2 and GM 2 /a*, respectively, where a* represents the radius of the 
PN spherical star defined by the conserved mass as 



a* = 




ao U + 



13GM* 
YlRc 2 



The normalized angular velocity and energy are written as 



Q 2 = 



2a2«3 

R 3 



1 + -^(2 -a\-a 



E = («2Q!3) 



1/3 



5R 2 
3 in 



9(a 2 a 3 ) 4/3 



2i? 4 



C s , 



(5-31) 



(5-32) 



1 1 

-~ + — ^(2 
2i? 10i? 3 



J? 2 



5Q!2Q!3 

55(a 2 Q 3 ) 2 / 3 



48i? 2 



5 a2a 3 
C s , 



4a 



+ «2 



(5-33) 



where -R = R/ci\ and C s = GM^j 'c 2 a*. 

Figures 15(a) and (b) show the energy as a function of the angular velocity in 
the corotating case (Jr = 0) and the irrotational one (fn = —2), respectively. The 
solid, dotted, short dashed, and long dashed lines represent the cases for C s = 0, 
0.01, 0.03 and 0.05, respectively. Note that for fn = —2, the minima correspond to 
the ISCO and for fn = 0, the minima correspond to the ISCCO. In Fig.16, we show 
fi at the ISCCO as a function of C s with (the dotted line) or without (the solid line) 
1PN effect of the self-gravity, where we use the normalization Q to present Fig.16 in 
order to compare it with Fig. 14. As expected the dotted line is above the solid line. 
Therefore from these calculations we can say that the 1PN effect of the self-gravity 
will increase fi and stabilize the binary. Note here that the Fig.16 corresponds to 
the lower left corner of Fig. 14. We checked the numerical value of Q in the solid 
lines of Fig. 14 and Fig.16 for the same relativity parameter and found the almost 
similar value. This is not strange. For small relativity parameter both treatments 
essentially agree since the pseudo-Newtonian potential is almost the same as the 
Newtonian potential. Therefore this result suggests that at least for small relativity 
parameter if we include 1PN effect of the self-gravity to TN, Q at ISCO (ISCCO) 
will increase. What will happen for large relativity parameter is not certain but we 
expect that the tendency will not change. 
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§6. Brief review on recent semi-relativistic works 

So far, we have investigated the location of the ISCO in the PN approximation. 
Although we have obtained qualitative properties of the GR correction for the lo- 
cation of the ISCO, the PN study is not sufficient to know the location accurately 
because it is accurate for small C B , but NS is a highly GR star of C s ~ 0.2. Thus, we 
need calculation including full GR terms or sufficient GR corrections. In the fully 
GR case, there is no equilibrium state of BNSs because of emission of GWs. Even 
so, we will be able to consider quasi-equilibrium states because tew is still much 
longer than the orbital period for orbits outside the ISCO. However, we have not 
known yet how to define the "quasi-equilibrium" in the fully GR case in contrast 
to the 1PN and 2PN cases in which we can distinguish quantities related to GWs 
from others. Thus, for fully GR calculation, we must begin from the unresolved first 
step, where we need to clarify the notion of the "quasi-equilibrium" and establish 
the formulation along the notion. 

Instead of fully GR one, two groups have been performing semi-GR calcula- 
tion. 1^ E§) The essence of their method is to assume that the three metric has the 
conformal flat form approximately. In this case, all the geometric variables which 
appear in the formalism are determined by solving Poisson equations, and wave equa- 
tions do not appear. Thus, effect of GWs is neglected and the formalism reduces to 
that similar to the PN one. Although some fully GR effects are taken into account 
in this formalism, the meaning of this approximation is still not sure. For example, 
if the orbital separation of BNSs is sufficiently large (i.e., tidal effects are not impor- 
tant) and spin angular velocity of each star is negligible, this formalism yields the 
exact GR solution. On the other hand, from the PN point of view, the formalism 
includes all the 1PN terms consistently, but not the 2PN terms when close BNSs are 
concerned. Thus, we may regard it as a formalism in which some GR corrections are 
merely included in the 1PN approximation. 

For a small compactness C s , this formalism will yield the same answer as that 
obtained by the 1PN calculation. Actually, this is the case for numerical results by 
Baumgarte et al. (BCSST) :© They obtained equilibrium states of corotating BNSs 
and their numerical results agree with those obtained in §4 for small C s cases, and 
even for large C s , they are similar to ours. From this fact, we may consider that 
this formalism is essentially very similar to the 1PN one. On the other hand, Wilson 
et al. did not obtain equilibrium states of corotating binary, but of non-uniform 
internal motion. Since the velocity field is different from the uniformly rotating one, 
the results may quantitatively disagree with ours and BCSST. However, qualitative 
nature also does not at all agree with ours and BCSST; in the analysis by Wilson et 
al., (l)i?isco decreases with increase of the GR effect and (2) the maximum density 
increases with decrease of the orbital separation, while the opposite relation holds 
in the analysis of ours and BCSST. There have been several researchers who doubt 
that accuxacv of their calculation is not enough and their results are something 
incorrect. 123 123113 However, we cannot completely deny a possibility that we miss 
some GR effects we have not been aware of. We may not reach a certain conclusion 
even for the case of this simplest version of semi-GR calculation at present. Hence, 
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to clarify the meaning of this approximation and accuracy of their numerical results, 
many works should be done in this area in near future. 

§7. Summary 

In this manuscript, we have presented recent results obtained from the study on 
the location of the ISCO in several levels of the PN approximation. The following is 
the brief summary. 

• Even in the Newtonian case, there exists the ISCO for binary of sufficiently stiff 
EOS. If the mass and the radius of each star are fixed, ^isco is larger for softer 
EOS. 

• There exist roughly two kinds of the 1PN corrections. One is the correction 
to the self-gravity of each star of binary, and the other is the correction to the 
gravity acting between two stars. The former one tends to increase i?iscO) but 
the latter one tends to decrease it. If we take into account both effects, however, 
the former effect is stronger than the latter one, and i?isco becomes large with 
increase of the 1PN correction. 

• The feature mentioned above is more remarkable for softer EOS when we fix 
the mass and radius. This is because for softer EOS, each star has the larger 
central density and is susceptible to the GR correction. 

• There has been no self consistent calculation including all the 2PN effects. But, 
there exist studies in which one only considers the 2PN gravity acting between 
two stars. In this case, L?isco is always smaller than that for the Newtonian 
case. If we include the PN effect of the self-gravity, L?isco will increase. 
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FIGURE CAPTIONS 



• Fig.l: Location of the ISCO is shown in each level of approximation, schemat- 
ically. 

• Fig. 2: Sketch of corotating binary. The origin of the coordinate is located at 
the center of mass of star 1. 

• Fig. 3: The energy(a) and the angular momentum(b) as functions of the angu- 
lar velocity for Newtonian incompressible binary stars of C / ' {GMao) 1 / 2 = —2, 
— 1.5, —1, —0.5, 0, 0.5, 1, 1.5, 2 and corotating binary. 



Fig. 4: The energy and the angular momentum as functions of the angular 
velocity for Newtonian corotating binary of compressible EOS. (a) for n = 0.25, 
(b) for n = 0.5, (c) for n = 0.75 and (d) for n = 1. 

Fig.5: J^N:iscco/^o as a function of the polytropic index n. Solid line denotes 
0. 248(1+0. 2n 1 - 5 ). Note that i^iscco was determined with ~ ±1% uncertainty. 
Fig. 6: The total energy(a) and angular momentum(b) of the equilibrium se- 
quence as functions of Q/Q* for 1PN incompressible binary stars. Solid, dot- 
ted, dashed and long dashed lines denote C s = GM*/c 2 a* = (the Newtonian 
case), 0.01, 0.03 and 0.05, respectively. 

Fig. 7: The orbital angular velocity at the ISCCO as a function of the compact- 
ness parameter C s = GM*/c 2 a* for 1PN incompressible binary stars. 
Fig. 8: The relation between R*/a* and J?/.!?* for n = 0.5 and C s = (open 
circles), 0.0327 (filled circles) and 0.0654 (filled squares). Dotted lines are drawn 
using Eq. ( |4-5| ) as counter parts of numerical results. 

Fig. 9: The energy and angular momentum as functions of QjQ* for n = 0.5(a), 
0.75(b) and 1(c) for the 1PN case. In Fig. 9(a), open circles, filled circles and 
open squares denotes sequences for C s = 0, 0.0327 and 0.0654, in Fig.9(b), 
they denote sequences for C s = 0, 0.02 and 0.04, and in Fig. 9(c), they denote 
sequences for C s = 0, 0.0167 and 0.0333, respectively. 

Fig. 10: The angular velocity at the ISCCO as a function of C s for n = 0.5(a), 
0.75(b) and 1(c). Note that ^iscco was determined with ~ ±1% uncertainty. 
Fig. 11: (a)The effective potentials of a test particle in the Schwarzschild metric 
(left) and the pseudo-Newtonian potential (right) as a function of the normal- 
ized distance R/r s . The vertical axes denote 

tf'Sch = yj(l ~ (2GM/Rc 2 ))((G 2 J| ch /i? 2 c 6 ) + 1) - 1 and ^ pscudo = -(GM/(R - 
r s )c 2 ) + (G 2 J 2 scudo /2R 2 c 6 ), respectively. The dots express the place of circular 
orbits. The value of the effective potential of the Schwarzschild black hole at 
the infinity is shifted to zero. (b)The fractional deviation of circular orbits of 
the pseudo-Newtonian potential from those of the Schwarzschild effective po- 
tential as a function of of the angular momentum. Circular orbits with different 
angular momentum related as J pS eudo = (9/8) 1 / 2 Js c h are compared (see text). 
Fig. 12: The total energy(a) and angular momentum(b) of the equilibrium se- 
quences as functions of the angular velocity in the case of corotating binaries. 
Solid, dotted, short dashed, long dashed, short dash-dotted, and long dash- 
dotted lines denote Gmi/aoc 2 = 0.05, 0.1, 0.15, 0.2 and 0.25, respectively. 
Fig. 13: The total energy(a) and angular momentum(b) for irrotational binaries 
as functions of the angular velocity. The conventions are the same as in Fig. 12. 
Fig. 14: The relation between the compactness of star 1 (Gmi/aoc 2 ) and the 
angular velocity of the binary Q = Q{Gm\/c') at the ISCO. Solid and dashed 
lines denote the cases of the corotating (Roche Ellipsoids) and irrotational bi- 
naries (Irrotational Roche-Riemann Ellipsoids), respectively. The value written 
at the upper right corner is that of the angular velocity in the point particle 
binary case. 
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• Fig. 15: The 1PN total energy without self-gravity terms as a function of the 
angular velocity in the corotating case(a) and the irrotational one(b), respec- 
tively. Solid, dotted, short dashed and long dashed lines represent the cases for 
C s = 0,0.01,0.03 and 0.05, respectively. 

• Fig.16: The relation between C s and Q = J?(GM*/c 3 ) at the ISCCO. Solid and 
dotted lines denote the cases of the 1PN corotating binaries without and with 
the 1PN self-gravity, respectively. 
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